miχpods: MOM6 dev#

  • daily cycle diagnostics

  • more wind dependencies figures + timing of descent marks

  • check background mixing formulation

  • bootstrap error on mean, median?

  • fix EUC max at 110

  • figure out time base

  • shear spectra

  • make χ, ε mean profiles.

  • do ε-Ri?

  • plot mean, medians on histograms.

Notes:#

Diagnostic notes/improvements#

  • match vertical resolution and extent

  • restrict TAO S2 to ADCP depth range

  • restrict to top 200m.

  • Filtering out MLD makes a big difference!

  • SInce we’re working with derivatives does the vertical grid matter (as long as it is coarser than observations)?

1 difvho ocean_vertical_heat_diffusivity 
2 difvso ocean_vertical_salt_diffusivity

Future model simulations#

  • Need lateral / neutral terms for total χ, ε

  • Add Atlantic TAO moorings + LOTUS, PAPA, others used in LMD94

  • Fix restart issue

Uncertainties#

  • Why can’t I reproduce the ε figure?

  • Why can’t I save the new station data

  • Need to consider Ri smoothing in models: From Gustavo:

    You mentioned some kind of inconsistency between your diagnostics and what MOM6 was doing for the interior shear mixing scheme. I am working to implement the option to apply vertical smoothing in Ri multiple times, instead of only once which is what we are doing right now. I noticed that the diagnostic ri_grad_shear is saved before the smoothing is applied. There is another diagnostic (ri_grad_shear_smooth) that saves the smoothed Ri. Perhaps you were looking at ri_grad_shear instead of ri_grad_shear_smooth and this can explain the inconsistency.

Diagnostics Effort notes#

  • go through prepare; pass in dataset and optional xgcm.Grid

  • different pre-processing steps for different datasets

  • order, sign of z coordinate is painful; here normalizing

  • composing with matplotlib subfigures

References#

Warner & Moum (2019)

Setup#

Hide code cell source 
%load_ext autoreload
%load_ext rich
%load_ext watermark

import warnings

import cf_xarray as cfxr
import dask
import dask_jobqueue
import dcpy
import distributed
import flox.xarray
import holoviews as hv
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import tqdm
import xgcm
import xrft
from datatree import DataTree
from mom6_tools.sections import preprocess_mom6_sections

import xarray as xr

%aimport pump
from pump import mixpods

hv.notebook_extension("bokeh")

cfxr.set_options(warn_on_missing_variables=False)
xr.set_options(keep_attrs=True, display_expand_data=False)

plt.style.use("bmh")
plt.rcParams["figure.dpi"] = 140

gcmdir = "/glade/campaign/cgd/oce/people/bachman/TPOS_1_20_20_year/OUTPUT/"  # MITgcm output directory
stationdirname = gcmdir

mixing_layer_depth_criteria = {
    "boundary_layer_depth": {"name": "KPPhbl|KPP_OBLdepth|ePBL_h_ML"},
}


from pump.catalog import catalog_dict

%watermark -iv
Hide code cell output 
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
The rich extension is already loaded. To reload it, use:
  %reload_ext rich
The watermark extension is already loaded. To reload it, use:
  %reload_ext watermark

dcpy         : 0.1.dev387+gd06c937
tqdm         : 4.65.0
re           : 2.2.1
pump         : 1.0+247.g1f1c5e1.dirty
xrft         : 0.0.0
xarray       : 2023.3.0
dask_jobqueue: 0.7.3
holoviews    : 1.15.4
intake       : 0.6.8
dask         : 2023.3.2
json         : 2.0.9
datatree     : 0.0.12
distributed  : 2023.3.2
flox         : 0.6.10
cf_xarray    : 0.8.0
numpy        : 1.23.5
matplotlib   : 3.7.1
ipywidgets   : 8.0.6
xgcm         : 0.6.1
sys          : 3.10.10 | packaged by conda-forge | (main, Mar 24 2023, 20:08:06) [GCC 11.3.0]
pandas       : 1.5.3
Hide code cell source 
if "client" in locals():
    client.close()
    del client
if "cluster" in locals():
    cluster.close()

cluster = dask_jobqueue.PBSCluster(
    cores=4,  # The number of cores you want
    memory="23GB",  # Amount of memory
    processes=1,  # How many processes
    queue="casper",  # The type of queue to utilize (/glade/u/apps/dav/opt/usr/bin/execcasper)
    local_directory="/local_scratch/pbs.$PBS_JOBID/dask/spill",
    log_directory="/glade/scratch/dcherian/dask/",
    resource_spec="select=1:ncpus=4:mem=23GB",  # Specify resources
    project="ncgd0011",  # Input your project ID here
    walltime="02:00:00",  # Amount of wall time
    interface="ib0",  # Interface to use
)
cluster.adapt(minimum_jobs=1, maximum_jobs=4)
client = distributed.Client(cluster)
client
Hide code cell output

Client

Client-51ce0365-ef4e-11ed-b816-3cecef1f3772

Connection method: Cluster object Cluster type: dask_jobqueue.PBSCluster
Dashboard: https://jupyterhub.hpc.ucar.edu/stable/user/dcherian/proxy/8787/status

Cluster Info

Read data#

LES#

%autoreload

les = mixpods.load_les_moorings()

--> The keys in the returned dictionary of datasets are constructed as follows:
	'latitude.longitude.month.kind.length'
100.00% [1/1 00:49<00:00]

microstructure#

micro = mixpods.load_microstructure()
micro

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:771: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:771: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:771: RuntimeWarning: divide by zero encountered in log10
  result_data = func(*input_data)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:771: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:771: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:771: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:771: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*

TAO#

tao_Ri = xr.load_dataarray(
    "tao-hourly-Ri-seasonal-percentiles.nc"
).cf.guess_coord_axis()

%autoreload

with dask.config.set(scheduler="threads"):
    tao_gridded = mixpods.load_tao()
tao_gridded

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:255: UserWarning: The specified Dask chunks separate the stored chunks along dimension "depth" starting at index 42. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:255: UserWarning: The specified Dask chunks separate the stored chunks along dimension "time" starting at index 199726. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:255: UserWarning: The specified Dask chunks separate the stored chunks along dimension "longitude" starting at index 2. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(

<xarray.Dataset>
Dimensions:             (time: 212302, depth: 61, depthchi: 6)
Coordinates: (12/19)
    deepest             (time) float64 dask.array<chunksize=(212302,), meta=np.ndarray>
  * depth               (depth) float64 -300.0 -295.0 -290.0 ... -10.0 -5.0 0.0
    eucmax              (time) float64 dask.array<chunksize=(33130,), meta=np.ndarray>
    latitude            float32 0.0
    longitude           float32 -140.0
    mld                 (time) float64 dask.array<chunksize=(212302,), meta=np.ndarray>
    ...                  ...
    oni                 (time) float32 -0.9 -0.9 -0.9 -0.9 ... nan nan nan nan
    en_mask             (time) bool False False False ... False False False
    ln_mask             (time) bool True True True True ... False False False
    warm_mask           (time) bool True True True True ... True True True True
    cool_mask           (time) bool False False False ... False False False
    enso_transition     (time) <U12 'La-Nina warm' ... '____________'
Data variables: (12/39)
    N2                  (time, depth) float64 dask.array<chunksize=(33130, 61), meta=np.ndarray>
    N2T                 (time, depth) float64 dask.array<chunksize=(33130, 61), meta=np.ndarray>
    Ri                  (time, depth) float64 dask.array<chunksize=(33130, 61), meta=np.ndarray>
    Rig_T               (time, depth) float64 dask.array<chunksize=(33130, 61), meta=np.ndarray>
    S                   (time, depth) float64 dask.array<chunksize=(33130, 61), meta=np.ndarray>
    S2                  (time, depth) float32 dask.array<chunksize=(33130, 61), meta=np.ndarray>
    ...                  ...
    Rig                 (time, depth) float64 dask.array<chunksize=(33130, 61), meta=np.ndarray>
    sst                 (time) float64 dask.array<chunksize=(33130,), meta=np.ndarray>
    Tflx_dia_diff       (time, depthchi) float64 nan nan nan nan ... nan nan nan
    ν                   (time, depthchi) float64 dask.array<chunksize=(33130, 6), meta=np.ndarray>
    Jb                  (time, depthchi) float64 dask.array<chunksize=(33130, 6), meta=np.ndarray>
    Rif                 (time, depthchi) float64 dask.array<chunksize=(33130, 6), meta=np.ndarray>
Attributes:
    CREATION_DATE:                23:26 24-FEB-2021
    Data_Source:                  Global Tropical Moored Buoy Array Project O...
    File_info:                    Contact: Dai.C.McClurg@noaa.gov
    Request_for_acknowledgement:  If you use these data in publications or pr...
    _FillValue:                   1.0000000409184788e+35
    array:                        TAO/TRITON
    missing_value:                1.0000000409184788e+35
    platform_code:                0n165e
    site_code:                    0n165e
    wmo_platform_code:            52321
np.log10(tao_gridded.eps).sel(time=slice("2005", "2015")).resample(
    time="M"
).mean().hvplot.quadmesh(clim=(-7.5, -6))

sub.v.hvplot()

sub = tao_gridded.sel(time="2010-01")

t = sub.theta.hvplot.quadmesh(cmap="turbo_r")
dt = (
    sub.theta - sub.theta.reset_coords(drop=True).cf.sel(Z=[0, -5]).cf.max("Z")
).hvplot.quadmesh(clim=(-0.15, 0.15), cmap="RdBu_r")
newmld = mixpods.get_mld_tao_theta(sub.theta.reset_coords(drop=True))

(
    dt
    * sub.reset_coords().mldT.hvplot.line(color="w", line_width=2)
    * newmld.reset_coords(drop=True).hvplot.line(color="orange", line_width=1)
).opts(frame_width=1200)

(
    tao_gridded.reset_coords().mldT.resample(time="5D").mean().hvplot.line()
    * mixpods.get_mld_tao_theta(tao_gridded.reset_coords().theta)
    .resample(time="5D")
    .mean()
    .hvplot.line()
)

tao_gridded.u.cf.plot()
tao_gridded.eucmax.plot()

[<matplotlib.lines.Line2D>]
_images/3a946f1b996d9cfc93e09718d8f29becf2633b045900e7ad003ba587cad536e4.png

MITgcm stations#

stations = pump.model.read_stations_20(stationdirname)

gcmeq = stations.sel(
    longitude=[-155.025, -140.025, -125.025, -110.025], method="nearest"
)
# enso = pump.obs.make_enso_mask()
# mitgcm["enso"] = enso.reindex(time=mitgcm.time.data, method="nearest")

# gcmeq["eucmax"] = pump.calc.get_euc_max(gcmeq.u)
# pump.calc.calc_reduced_shear(gcmeq)

# oni = pump.obs.process_oni()
# gcmeq["enso_transition"] = mixpods.make_enso_transition_mask(oni).reindex(time=gcmeq.time.data, method="nearest")


mitgcm = gcmeq.sel(longitude=-140.025, method="nearest")
metrics = mixpods.normalize_z(pump.model.read_metrics(stationdirname), sort=True)
mitgcm = mixpods.normalize_z(mitgcm, sort=True)

mitgcm_grid = xgcm.Grid(
    metrics.sel(latitude=mitgcm.latitude, longitude=mitgcm.longitude, method="nearest"),
    coords=({"Z": {"center": "depth", "outer": "RF"}, "Y": {"center": "latitude"}}),
    metrics={"Z": ("drF", "drC")},
    periodic=False,
    boundary="fill",
    fill_value=np.nan,
)

mitgcm.theta.attrs["standard_name"] = "sea_water_potential_temperature"
mitgcm.salt.attrs["standard_name"] = "sea_water_salinity"
mitgcm["KPPviscAz"].attrs["standard_name"] = "ocean_vertical_viscosity"
mitgcm["KPPdiffKzT"].attrs["standard_name"] = "ocean_vertical_heat_diffusivity"

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/indexing.py:1374: PerformanceWarning: Slicing is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array[indexer]

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array[indexer]
  return self.array[key]
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/indexing.py:1374: PerformanceWarning: Slicing is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array[indexer]

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array[indexer]
  return self.array[key]

mitgcm = (
    mixpods.prepare(mitgcm, grid=mitgcm_grid, oni=pump.obs.process_oni())
    .sel(latitude=0, method="nearest")
    .cf.sel(Z=slice(-250, 0))
)

mitgcm_grid

<xgcm.Grid>
Z Axis (not periodic, boundary='fill'):
  * center   depth --> outer
  * outer    RF --> center
Y Axis (not periodic, boundary='fill'):
  * center   latitude

mitgcm.u.cf.plot()
mitgcm.mldT.reset_coords(drop=True).cf.plot()
mitgcm.eucmax.reset_coords(drop=True).cf.plot()

[<matplotlib.lines.Line2D>]
_images/d46c7b8c2987c0ef2648ef1244623ce16ba4db3433f18083967672403b8ba06c.png 
mixpods.plot_n2s2pdf(mitgcm.n2s2pdf.sel(enso_transition_phase="none"))

<matplotlib.contour.QuadContourSet>
_images/0687c893e1eaca992b13f8fe6166eda6df0e63d890df04891f6f05f5cee52932.png

MOM6#

Generate kechunk JSONs#

catalog_sub = {k: v for k, v in catalog_dict.items() if k}
catalog_sub

{
    'baseline.001': ('baseline', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods'),
    'epbl': ('ePBL', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods'),
    'kpp.lmd.002': (
        'KPP Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.002.mixpods'
    ),
    'kpp.lmd.003': (
        'KPP Ri0=0.5, Ric=0.2,',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.003.mixpods'
    ),
    'kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods'
    ),
    'baseline.N150': ('baseline N=150', 'gmom.e23.GJRAv3.TL319_t061_zstar_N150.baseline.mixpods'),
    'kpp.lmd.004.N150': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5, N=150',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N150.kpp.lmd.004.mixpods'
    ),
    'baseline.hb': ('baseline', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.hb'),
    'new_baseline.hb': ('KD=0, KV=0', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb'),
    'new_baseline.kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.004.mixpods'
    ),
    'new_baseline.kpp.lmd.005': (
        'KPP ν0=2.5, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods'
    )
}

%autoreload

from mom6_tools.kerchunk import generate_references_for_stream

for _, casename in tqdm.tqdm(catalog_sub.values()):
    caseroot = f"{mixpods.ROOT}/cesm/{casename}"
    if "N150" in casename:
        continue
    print(caseroot)
    for stream in ["h", "hm", "hm.wci", "sfc"]:
        generate_references_for_stream(
            caseroot=caseroot,
            stream=stream,
            missing_stream="warn",
            existing_output="overwrite",
        )

../pump-catalog catalog with 4 dataset(s) from 4 asset(s):

unique
casename 1
stream 4
path 4
baseline 1
levels 1
frequency 1
variables 75
shortname 1
description 1
derived_variables 0 
  0%|          | 0/1 [00:00<?, ?it/s]

/glade/campaign/cgd/oce/projects/pump//cesm/gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods

/glade/u/home/dcherian/python/mom6-tools/mom6_tools/kerchunk.py:104: RuntimeWarning: No files found for caseroot: /glade/campaign/cgd/oce/projects/pump//cesm/gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods, stream: hm.wci
  warnings.warn(f"No files found for caseroot: {caseroot}, stream: {stream}", RuntimeWarning)
100%|██████████| 1/1 [01:34<00:00, 94.96s/it]

import pathlib

root = pathlib.Path(
    "/glade/campaign/cgd/oce/projects/pump//cesm/gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods"
)

%autoreload

from mom6_tools.kerchunk import combine_stream_jsons_as_groups

for casename in tqdm.tqdm(catalog_sub.df["casename"].unique()):
    caseroot = f"{mixpods.ROOT}/cesm/{casename}"
    combine_stream_jsons_as_groups(caseroot=caseroot, streams=None)

100%|██████████| 10/10 [00:13<00:00,  1.37s/it]

staticfile = (
    "/glade/u/home/dcherian/campaign-oce/projects/pump/cesm/gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb/"
    "run/gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb.mom6.static.nc"
)

import fsspec

from mom6_tools.kerchunk import open_references_as_xarray

fs = fsspec.filesystem(
    "reference",
    fo=f"{caseroot}/run/jsons/combined.json",
    skip_instance_cache=True,
)
mapper = fs.get_mapper(root="")

import datatree

xr.open_dataset(mapper, engine="zarr", use_cftime=True, consolidated=False)

<xarray.Dataset>
Dimensions:  ()
Data variables:
    *empty*
open_references_as_xarray(f"{caseroot}/run/jsons/combined.json")

<xarray.Dataset>
Dimensions:  ()
Data variables:
    *empty*

calculate ONI#

(
    oniobs.hvplot.line(x="time", label="obs")
    * onimom6.hvplot.line(x="time", label="MOM6")
).opts(ylabel="ONI [°C]")

oniobs.enso_transition_mask.plot()
onimom6.enso_transition_mask.plot(color="r")

[<matplotlib.lines.Line2D>]
_images/94103287fcd3b60c9257f03a6985bb4b98cad673ba31d02f429aca0a3f467360.png

MOM6 sections#

Combine sections#
catalog

{
    'baseline': ('baseline', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods'),
    'epbl': ('ePBL', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods'),
    'kpp.lmd.002': (
        'KPP Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.002.mixpods'
    ),
    'kpp.lmd.003': (
        'KPP Ri0=0.5, Ric=0.2,',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.003.mixpods'
    ),
    'kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods'
    ),
    'baseline.N150': ('baseline N=150', 'gmom.e23.GJRAv3.TL319_t061_zstar_N150.baseline.mixpods'),
    'kpp.lmd.004.N150': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5, N=150',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N150.kpp.lmd.004.mixpods'
    ),
    'new_baseline.hb': ('KD=0, KV=0', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb'),
    'new_baseline.kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.004.mixpods'
    ),
    'new_baseline.kpp.lmd.005': (
        'KPP ν0=2.5, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods'
    )
}

casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods"
ds = mixpods.read_mom6_sections(casename)

100%|██████████| 39/39 [02:32<00:00,  3.90s/it]

ds.drop_vars(
    ["average_DT", "average_T2", "average_T1", "time_bnds"], errors="ignore"
).chunk({"time": 24 * 365}).to_zarr(
    f"/glade/scratch/dcherian/archive/{casename}/ocn/moorings/tao.zarr",
    consolidated=True,
    mode="w",
)

<xarray.backends.zarr.ZarrStore object at 0x2aafa7cb0c10>

reload = mixpods.load_mom6_sections(casename)
reload.uo.count().compute()

%autoreload

casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods"
casename = mixpods.mom6_sections_to_zarr(casename)

100%|██████████| 33/33 [00:55<00:00,  1.69s/it]

%autoreload

casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.004.mixpods"
casename = mixpods.mom6_sections_to_zarr(casename)

100%|██████████| 33/33 [00:59<00:00,  1.81s/it]

casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N150.baseline.mixpods"
mixpods.mom6_sections_to_zarr(casename)

100%|██████████| 39/39 [00:29<00:00,  1.30it/s]

casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods"
mixpods.mom6_sections_to_zarr(casename)

casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods"
mixpods.mom6_sections_to_zarr(casename)

100%|██████████| 31/31 [01:50<00:00,  3.57s/it]
/glade/u/home/dcherian/pump/pump/mixpods.py:887: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar.  This may lead to subtle errors in operations that depend on the length of time between dates.
  mom6tao["time"] = mom6tao.indexes["time"].to_datetimeindex()

casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.002.mixpods"
mixpods.mom6_sections_to_zarr(casename)

casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.003.mixpods"
mixpods.mom6_sections_to_zarr(casename)

100%|██████████| 38/38 [02:07<00:00,  3.36s/it]
/glade/u/home/dcherian/pump/pump/mixpods.py:900: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar.  This may lead to subtle errors in operations that depend on the length of time between dates.
  mom6tao["time"] = mom6tao.indexes["time"].to_datetimeindex()

casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods"
mixpods.mom6_sections_to_zarr(casename)

Read sections#

dask.config.set(scheduler=client)
m = mixpods.read_mom6_sections(
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods"
)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:274: FutureWarning: env_extra has been renamed to job_script_prologue. You are still using it (even if only set to []; please also check config files). If you did not set job_script_prologue yet, env_extra will be respected for now, but it will be removed in a future release. If you already set job_script_prologue, env_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/pbs.py:82: FutureWarning: project has been renamed to account as this kwarg was used wit -A option. You are still using it (please also check config files). If you did not set account yet, project will be respected for now, but it will be removed in a future release. If you already set account, project is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
100%|██████████| 21/21 [02:14<00:00,  6.42s/it]
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/coding/times.py:360: FutureWarning: Index.ravel returning ndarray is deprecated; in a future version this will return a view on self.
  sample = dates.ravel()[0]
/glade/u/home/dcherian/pump/pump/mixpods.py:847: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar.  This may lead to subtle errors in operations that depend on the length of time between dates.
  mom6tao["time"] = mom6tao.indexes["time"].to_datetimeindex()
/glade/u/home/dcherian/pump/pump/mixpods.py:859: UserWarning: Kv_v not present. Assuming equal to Kv_u
  warnings.warn("Kv_v not present. Assuming equal to Kv_u")

m.drop_vars(["average_DT", "average_T1", "average_T2", "time_bnds"]).chunk(
    {"time": 365 * 24}
).to_zarr(
    "/glade/scratch/dcherian/archive/gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods/ocn/moorings/tao.zarr",
    mode="w",
)

m.sel(yh=0, method="nearest")[["ePBL_h_ML", "mlotst"]].to_array().hvplot.line(
    by="variable", x="time"
)

m.sel(yh=0, method="nearest")[["Kd_heat", "Kd_ePBL"]].to_array().hvplot.line(
    by="variable", groupby="time", logy=True, ylim=(1e-6, 1e-1), xlim=(0, 500)
)

mom6140 = mixpods.load_mom6_sections(casename).load()

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1471: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1471: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1471: UserWarning: Metric at ('time', 'zi', 'yq', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1638: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1471: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1471: UserWarning: Metric at ('time', 'zi', 'yq', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(

mom6140.uo.cf.plot(robust=True)
mom6140.eucmax.cf.plot()
mom6140.mldT.cf.plot(lw=0.5, color="k")

mixpods.plot_n2s2pdf(mom6140.n2s2pdf.sel(enso_transition_phase="none"))

<matplotlib.contour.QuadContourSet>
_images/b8bed40b81209ef9ccde157be74e5b3a06571f6e7cfba20a8759504585d48dfe.png

Pick simulations#

Build dict#

import warnings

datasets = {
    # "TAO": tao_gridded,
    # "MITgcm": mitgcm,
}

catalog_sub = {
    k: catalog_dict[k]
    for k in catalog_dict.keys()
    if (k == "kpp.lmd.004") or ("baseline" in k and "150" not in k)
}
display(catalog_sub)

{
    'baseline': ('baseline', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods'),
    'kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods'
    ),
    'new_baseline.hb': ('KD=0, KV=0', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb'),
    'new_baseline.kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.004.mixpods'
    ),
    'new_baseline.kpp.lmd.005': (
        'KPP ν0=2.5, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods'
    )
}

%autoreload

with warnings.catch_warnings():
    warnings.simplefilter("ignore", category=UserWarning)
    warnings.simplefilter("ignore", category=FutureWarning)
    for short_name, (long_name, folder) in tqdm.tqdm(catalog_sub.items()):
        datasets.update(
            {
                short_name: mixpods.load_mom6_sections(folder).assign_attrs(
                    {"title": long_name}
                )
            }
        )

100%|██████████| 5/5 [00:35<00:00,  7.16s/it]

datasets["TAO"] = DataTree(tao_gridded)

datasets["les"] = les["0.0.-140.oct.average.month"].to_dataset()

# Offset LES to work with slicing below
datasets["les"]["time"] = datasets["les"]["time"] + pd.Timedelta(days=25 * 365)

Build tree#

tree = DataTree.from_dict(datasets)

tree = tree.sel(time=slice("2000", "2017"))

ref = tree["TAO"].ds.reset_coords(drop=True).cf.sel(Z=slice(-120, None))
counts = np.minimum(ref["S2"].cf.count("Z"), ref["N2T"].cf.count("Z")).load()


def calc_histograms(ds):
    ds = ds.copy()
    ds["tao_mask"] = counts.reindex(time=ds.time, method="nearest") > 5
    ds["tao_mask"].attrs = {
        "description": "True when there are more than 5 5-m T, u, v in TAO dataset"
    }
    # ds = ds.where(ds.tao_mask)
    return ds.update(mixpods.pdf_N2S2(ds))


tree = tree.map_over_subtree(calc_histograms)
if "les" in tree:
    tree["les"] = tree["les"].isel(z=slice(-2))
    tree["les"]["KT"].attrs["standard_name"] = "ocean_vertical_heat_diffusivity"

if "micro" in locals():
    tree.update(micro)

tree

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*
%autoreload

dailies = tree.map_over_subtree(mixpods.daily_composites)
dailies

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out,

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*
tree = mixpods.persist_tree(tree)

mixpods.validate_tree(tree)

Comparing shear prod to eps#

(ds.shear_prod.mean("time").hvplot.line() * ds.eps.mean("time").hvplot.line()).opts(
    hv.opts.Curve(invert_axes=True, ylim=(1e-10, None), logx=True)
)

Daily composites#

  • S2-N2 histogram by hour of day?

  • daily composites by el-nino phase

  • figure out what wind stress is used

    • wind stress diurnal cycle

    • wind stress in model is relative wind

    • look at the coupled model

  • either turbulence parameters

    • Jess Masich’s papers?

  • MOM6 uses Large & Pond; χpods use Large & Pond

    • does it have gustiness; MOM6 uses wind averaged over coupling interval;

  • Look at 69m-89m range with medium wind stress

  • Need to get rid of EUC effect:

    • either EUC-relative coordinate system

    • or use K, which removes a S2 factor

  • Do buoyancy flux or Γ

ε, τ#

  • How well does the model get τ

    • τ bins are [0, 0.04, 0.075, inf]

def groupby_mean(node):
    return node["tau"].groupby_bins(node.tau, bins=np.linspace(0, 0.2, 51)).count()


tau_hist = (
    tree.dc.subset_nodes(["tau"])
    .map_over_subtree(groupby_mean)
    .reset_coords(drop=True)
    .dc.concatenate_nodes()
)
tau_hist["tau_bins"] = pd.IntervalIndex(tau_hist.tau_bins.data).mid

tau_hist.hvplot.step(by="node")

bins = np.linspace(0, 0.2, 51)
for name, node in tree.children.items():
    node["tau"].plot.hist(bins=bins, label=name, histtype="step", lw=3, density=True)
plt.legend()
dcpy.plots.linex([0.04, 0.075])
_images/966e87b1b06744b5507cb79d782543af57dcd110a6d571c5a1db1ab4c1e5261a.png

Verify depth is normalized#

for name, ds in datasets.items():
    (ds.cf["sea_water_x_velocity"].cf["Z"].plot(marker=".", ls="none", label=name))
plt.legend()

<matplotlib.legend.Legend>
_images/2372712bb196933da3d938a29d3f04e0e217fe29d04bb036dbceb08ab9ce12da.png

Compare EUC maximum and MLD#

Monthly climatology#

import tqdm

for node in tqdm.tqdm(tree.children):
    tree[node]["mldT"] = tree[node]["mldT"].load()

100%|██████████| 3/3 [00:28<00:00,  9.35s/it]

def to_dataset(tree):
    concat = xr.concat([v.ds for v in tree.children.values()], dim="node")
    concat["node"] = list(tree.children.keys())
    return concat


clim = to_dataset(
    tree.map_over_subtree(
        lambda x: x.reset_coords()[["mldT", "eucmax"]].groupby("time.month").mean()
    )
).load()
clim

<xarray.Dataset>
Dimensions:  (node: 3, month: 12)
Coordinates:
  * month    (month) int64 1 2 3 4 5 6 7 8 9 10 11 12
  * node     (node) <U16 'TAO' 'kpp.lmd.004' 'kpp.lmd.004.N150'
Data variables:
    mldT     (node, month) float64 -24.54 -16.38 -12.75 ... -21.87 -25.93 -31.85
    eucmax   (node, month) float64 -107.1 -100.9 -93.21 ... -103.6 -104.7 -103.6
hv.Layout(
    [
        hv.Overlay(
            [
                tree[node]
                .ds["mldT"]
                .reset_coords(drop=True)
                .groupby("time.month")[month]
                .hvplot.hist(label=node, legend=True)
                .opts(frame_width=150)
                for node in tree.children
            ]
        ).opts(title=str(month))
        for month in np.arange(1, 13)
    ]
).cols(4)

(clim.mldT.hvplot(by="node") + clim.eucmax.hvplot(by="node")).cols(1)

mixpods.plot_timeseries(tree, "mldT", obs="TAO")

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/bokeh/core/property/bases.py:259: DeprecationWarning: elementwise comparison failed; this will raise an error in the future.
  return new == old

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/bokeh/core/property/bases.py:259: DeprecationWarning: elementwise comparison failed; this will raise an error in the future.
  return new == old

datasets["TAO"].eucmax.attrs["long_name"] = "EUC maximum"
mixpods.plot_timeseries(tree.sel(time="2008"), "eucmax", obs="TAO")

MLD Maps#

def read_sfc(casename):
    path = f"/glade/scratch/dcherian/archive/{casename}/ocn/hist/*sfc*00[4-7]*.nc"
    sfc = xr.open_mfdataset(
        path, data_vars="minimal", coords="minimal", compat="override", parallel=True
    )
    return sfc


sfc = DataTree()
casenames = {
    "MOM6 KPP": "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods",
    "MOM6 ePBL": "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods",
}

for name, casename in casenames.items():
    sfc[name] = DataTree(read_sfc(casename))

sfc

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*
mldclim = sfc.map_over_subtree(
    lambda node: node.mlotst.groupby("time.month").mean()
).compute()

filepath = (
    "/glade/work/gmarques/cesm/datasets/MLD/deBoyer/deBoyer_MLD_remapped_to_tx06v1.nc"
)
mldclim["obs"] = DataTree(xr.open_dataset(filepath))
mldclim

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/coding/times.py:699: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/indexing.py:524: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(array[self.key], dtype=None)

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*
mldclim.

mldclim["MOM6 KPP"] = DataTree(
    read_sfc("gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods")
    .mlotst.groupby("time.month")
    .mean()
)

mldclim["MOM6 ePBL"] = DataTree(
    read_sfc("gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods")
    .mlotst.groupby("time.month")
    .mean()
)

(
    (mlotst - mldclim.mld.groupby("time.month").first())
    .sel(yh=slice(-5, 5), xh=slice(-170, -80))
    .plot(col="month", robust=True, col_wrap=3, aspect=2)
)

<xarray.plot.facetgrid.FacetGrid>
_images/379b40bc88e48b23199166fe6ba7e8a7925e86dce6842cc7315c2a11db8e079a.png

Mixing layer depth?#

todrop = ["TAO"]
notao = DataTree.from_dict(
    {k: v.ds for k, v in tree.children.items() if k not in todrop}
)
with cfxr.set_options(custom_criteria=mixing_layer_criteria):
    plot = mixpods.plot_timeseries(notao, "boundary_layer_depth", obs=None)
plot

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(

Mean profiles: mean, std#

Model mean S2 is lower by a factor of 2

S2 = mixpods.plot_profile_fill(tree, "S2", "S^2")
N2 = mixpods.plot_profile_fill(tree, "N2T", "N_T^2")
u = mixpods.plot_profile_fill(tree, "sea_water_x_velocity", "u")
T = mixpods.plot_profile_fill(tree, "sea_water_potential_temperature", "T")
Ri = mixpods.plot_median_Ri(tree)
Hide code cell source 
S = mixpods.plot_profile_fill(tree, "sea_water_salinity", "S")
v = mixpods.plot_profile_fill(tree, "sea_water_y_velocity", "v")
# dTdz = plot_profile_fill(tree, "dTdz", "∂T/∂z")
# dSdz = plot_profile_fill(tree, "dSdz", "∂S/∂z")

((S2 + N2).opts("Curve", xrotation=30) + Ri.opts(xaxis="top", xrotation=30)).opts(
    sizing_mode="stretch_width"
)

((S2 + N2).opts("Curve", xrotation=30) + Ri.opts(xaxis="top", xrotation=30)).opts(
    sizing_mode="stretch_width"
)

(T + u.opts(ylim=(-0.5, 1.5))).opts(sizing_mode="stretch_width")

(T + S + u.opts(ylim=(-0.5, 1.5)) + v.opts(ylim=(-0.5, 0.5))).opts(
    sizing_mode="stretch_width"
)

S2 = plot_profile_fill(tree, "S2", "S^2")
u = plot_profile_fill(tree, "sea_water_x_velocity", "u")

popts = (
    hv.opts.Curve(fontscale=1.5, line_width=3, color=hv.Cycle("Dark2")),
    hv.opts.Area(fill_color=hv.Cycle("Dark2")),
)

h = (S2.opts(show_legend=False) + u).opts(
    hv.opts.Curve("Curve", xaxis="bottom", xlim=(-200, 0)),
    hv.opts.Overlay(frame_width=400, frame_height=600),
    *popts,
)
h

Remap to EUC coordinate#

gcmeq.coords["zeuc"] = gcmeq.depth - gcmeq.eucmax
remapped = flox.xarray.xarray_reduce(
    gcmeq.drop_vars(["SSH", "KPPhbl", "mld", "eucmax"], errors="ignore"),
    "zeuc",
    dim="depth",
    func="mean",
    expected_groups=(np.arange(-250, 250.1, 5),),
    isbin=(True,),
)
remapped["zeuc_bins"].attrs["units"] = "m"
remapped

<xarray.Dataset>
Dimensions:          (time: 174000, longitude: 4, zeuc_bins: 100)
Coordinates:
    latitude         float64 0.025
  * longitude        (longitude) float64 -155.0 -140.0 -125.0 -110.0
  * time             (time) datetime64[ns] 1998-12-31T18:00:00 ... 2018-11-06...
    cool_mask        (time) bool True True True True ... False False False False
    warm_mask        (time) bool False False False False ... True True True True
  * zeuc_bins        (zeuc_bins) object (-250.0, -245.0] ... (245.0, 250.0]
Data variables: (12/25)
    DFrI_TH          (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    KPPdiffKzT       (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    KPPg_TH          (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    KPPviscAz        (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    VISrI_Um         (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    VISrI_Vm         (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    ...               ...
    S2               (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    shear            (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    N2               (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    shred2           (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    Ri               (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    enso_transition  (zeuc_bins, time) <U12 'La-Nina cool' ... 'El-Nino warm'
Attributes:
    easting:   longitude
    northing:  latitude
    title:     Station profile, index (i,j)=(1201,240)
remapped = remapped.persist()

cluster.scale(6)

Seasonal median Ri: (0, 140)#

remapped["longitude"] = [-155.0, -140, -125, -110]

fg = (
    remapped.Ri.groupby("time.season")
    .median()
    .sel(season=["DJF", "MAM", "JJA", "SON"], longitude=[-140, -110])
    .cf.plot(
        col="season", row="longitude", xlim=(0, 1.5), ylim=(-20, None), label="MITgcm"
    )
)
fg.data = tao_Ri.cf.sel(quantile=0.5, longitude=[-140, -110])
fg.map_dataarray_line(
    xr.plot.line, x=None, y="zeuc", hue="season", color="k", label="TAO"
)
fg.map(lambda: plt.axvline(0.25))
fg.axes[-1, -1].legend(bbox_to_anchor=(1.5, 1))

<matplotlib.legend.Legend>
_images/de806fb2b38296bab1ec713f0632b54dff28a3fa7f59f43ee018acdefada8e19.png

1D Climatological heat budget#

  • Using TAO swnet, qnet measurements instead of Tropflux

tree["TAO"] = DataTree(tao_gridded)

mixpods.climo_heat_budget_1d(tree["TAO"].ds, penetration="moum")
mixpods.climo_heat_budget_1d(tree["TAO"].ds, penetration="mom")

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
_images/c81627abe7b7b8ad053498ae040d320f7f62ea4d9a4c7708bc775e490af071ae.png _images/8d3bf9bad20c6663b5023fd3acb30fbbbdc5c71561ce3d0818e2747e0166868d.png 
f, ax = plt.subplots(
    1,
    len(tree),
    constrained_layout=True,
    squeeze=False,
    sharex=True,
    sharey=True,
    figsize=(10, 3),
)

for axx, (name, node) in zip(ax.flat, tree.children.items()):
    mixpods.plot_climo_heat_budget_1d(node.ds, mxldepth=-40, penetration="mom", ax=axx)
    axx.set_title(name)

dcpy.plots.clean_axes(ax)

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
_images/fb5f66b090672708c55e4945761e9ee7c7f1d7203227dd8464c79bac381e3988.png 
mixpods.plot_climo_heat_budget_1d(tree["kpp.lmd.004.N150"].ds, penetration="mom")

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
_images/bb12a5ffc0d1660adedecbe434c2486c1dd78b6d4d48551b2d0ae7ac2c818880.png

Stability diagram: 4N2-S2 plane#

Warner & Moum (2019):

Both velocity and temperature are hourly averaged before interpolation to 5-m vertical bins

Contours enclose 50% of data

Compare daily vs hourly frequency#

resampled = tree.from_dict(
    {k: v.ds.resample(time="D").mean() for k, v in tree.children.items()}
)
# resampled["TAO"].update(mixpods.prepare(resampled["TAO"].ds, oni=mixpods.process_oni()))
for name, node in resampled.children.items():
    node.update(
        mixpods.prepare(
            node.ds,
            oni=tree[name]["oni"].resample(time="D").mean().reset_coords(drop=True),
        )
    )
    node.update(mixpods.pdf_N2S2(node.ds))
    mixpods.validate_tree(node)

resampled = mixpods.load_tree(resampled)

tree = mixpods.load_tree(tree)

fig = plt.figure(constrained_layout=True, figsize=(12, 6))

subfigs = fig.subfigures(2, 1)
mixpods.plot_stability_diagram_by_dataset(tree, fig=subfigs[0])
mixpods.plot_stability_diagram_by_dataset(resampled, fig=subfigs[1])
_images/e3c72d6db28d6cceaeca1332dcdea39cd873540b9b68965123d11396ffb7ef33.png

Older#

Hide code cell source 
mixpods.plot_stability_diagram_by_dataset(datasets)
_images/2027c31463a49bfd091b264ec6d31d57a6debcc9415b3587f08d751aaac9ded3.png
Hide code cell source 
mixpods.plot_stability_diagram_by_phase(datasets, obs="TAO", fig=None)
_images/311676930e479f61019987704cb52d24a341c248f3c24ad42f0aacf6347cda46.png

Composed#

excluding MLD#

  • have not matched vertical grid spacing yet.

  • minor axis is smaller for increasing shear with KPP

  • prandtl number

  • interpolate to TAO grid

  • hybrid HYCOM-like case may have equatorial refinement of vertical grid above EUC

  • compare hybrid coordinate case resolution

  • instrumental error on vertical gradients for S2, N2;

    • if added to model, would it change the bottom tail.

  • Horizontal or vertical res:

    • MOM6 at 5km? hah

    • TIW compositing? La-NIna dec?

      • top end’s are similar, so maybe not?

  • sensitivity of cold tongue mean state

    • horizontal viscosity; controls EUC strength? (increase?)

    • vertical res; yan li’s paper

    • POP :thermocline always too diffuse and EUC too thick;

    • higher order interpolants can control amount diffusion when remapping

  • get S2, N2 from TPOSE?

Questions:

  1. internal wave spectrum?

  2. minor axis is smaller for model, particularly MITgcm-> too diffusive?

  3. Can we think of this as the model relaxing to critical Ri too quickly

tree = mixpods.load_tree(tree)

fig = plt.figure(constrained_layout=True, figsize=(12, 6))

subfigs = fig.subfigures(2, 1)
mixpods.plot_stability_diagram_by_dataset(tree, fig=subfigs[0])
mixpods.plot_stability_diagram_by_phase(tree, fig=subfigs[1])
Hide code cell source 
fig = plt.figure(constrained_layout=True, figsize=(12, 6))

subfigs = fig.subfigures(2, 1)
mixpods.plot_stability_diagram_by_dataset(tree, fig=subfigs[0])
mixpods.plot_stability_diagram_by_phase(tree, fig=subfigs[1])
_images/050e916a870137e4d90530d173b39cbc75d84ebadff2f5484a33570e58c658a2.png

including MLD#

Hide code cell source 
fig = plt.figure(constrained_layout=True, figsize=(12, 6))
subfigs = fig.subfigures(2, 1, height_ratios=[1.2, 1])
top = subfigs[0].subfigures(1, 3, width_ratios=[1, 5, 1])

mixpods.plot_stability_diagram_by_dataset(datasets, fig=top[1])
mixpods.plot_stability_diagram_by_phase(datasets, fig=subfigs[1])
_images/678886ada6c73c9a3c0baadcc9597ed60922d04ef61fb87e2dde1c0d8ec1bafa.png 
(
    tao_gridded.n2s2pdf.sel(
        enso_transition_phase=[
            "La-Nina cool",
            "La-Nina warm",
            "El-Nino warm",
            "El-Nino cool",
        ]
    )
    .sum("N2_bins")
    .plot(hue="enso_transition_phase")
)
plt.figure()
(
    tao_gridded.n2s2pdf.sel(
        enso_transition_phase=[
            "La-Nina cool",
            "La-Nina warm",
            "El-Nino warm",
            "El-Nino cool",
        ]
    )
    .sum("S2_bins")
    .plot(hue="enso_transition_phase")
)

[<matplotlib.lines.Line2D>,
 <matplotlib.lines.Line2D>,
 <matplotlib.lines.Line2D>,
 <matplotlib.lines.Line2D>]
_images/4e7b99af96d4ce455831441e3352c755dd0fed619ae4d068afa92e5d99bb25d8.png _images/1db17f0d2ed4f555c9f2e749345f58130bb2819db51a3f0a1ee17a981b967f7e.png 
oni = pump.obs.process_oni().sel(time=slice("2005-Oct", "2015"))
(
    oni.hvplot.step(color="k")
    + pump.obs.make_enso_transition_mask()
    .sel(time=slice("2005-Jun", "2015"))
    .reset_coords(drop=True)
    .hvplot.step()
).cols(1)

Turbulence#

Dan’s suggestions#

  1. I liked your effort to look at shear production= KmdU/dz dot dU/dz = dissipation. I mentioned that it might be worth looking at buoyancy flux = KT dT/dz * db/dT + KS dS/dz * db/dS with db/dT, db/dS from equation of state. I think I said KT*N^2, but I hope after checking that this is the same in this region (K_T~K_S). You could then divide buoyancy flux by 0.2 to estimate epsilon for comparison with your estimate based on shear production.

  1. Dimensionless parameters. We alluded to this, but I think it would be interesting to explicitly quantify flux Richardson number Rif= buoyancy flux/shear production (as defined above). And Prt=KT/Km=Rif/Rig, where Rig is the usual gradient Richardson number. Following a similar logic, you might look at Reb= epsilon/nu N^2 with nu constant 1e-6 molecular viscosity. We discussed looking at diffusivity and viscosity in addition to epsilon. But note the similarity between dimensionless diffusivity KT/kappa and viscosity Km/nu, i.e. Reb=Rig*KM/nu. (with kappa a constant molecular diffusivity to make dimensionless)

  1. I like how you continue to push what we can learn from the distributions in S2-N2 parameter space in addition to physical space. It is good for comparing with limited observations or process models. I think you can push this even further to develop metrics quantifying how the momentum and buoyancy fluxes “tango” with each other and with S2 and N2. In addition, perhaps it would be valuable to flip the perspective and look at shear-production-buoyancy flux space? Or consider various two-dimensional spaces of non-dimensional parameters mentioned above: Rig, Reb, Rif, KT/kappa, Km/nu.

  1. Finally, I think it would be interesting to try to examine how momentum flux and shear production below the mixed layer relate to wind stress as well as the local shear. There might be a relationship that emerges in KPP, even if the K profile is formally dependent only on local Rig at these depths.

Daily composites#

h = mixpods.plot_daily_composites(
    dailies.dc.subset_nodes(["S2", "N2", "Rig_T"]), logy=False
).opts(
    hv.opts.Curve(frame_width=200, xrotation=45),
    hv.opts.GridSpace(show_legend=True),
)
h

h = mixpods.plot_daily_composites(
    dailies.dc.subset_nodes(["eps", "KT", "chi"]), logy=True
)
h.opts(
    hv.opts.Curve(frame_width=200, xrotation=45),
    hv.opts.GridSpace(show_legend=True),
)

ε-Ri means#

mixpods.map_hvplot(
    lambda dt, name, muted: mixpods.plot_eps_ri_hist(
        dt["eps_ri"].load(), label=name, muted=muted
    ),
    tree.children,
).opts(show_legend=True).opts(hv.opts.Curve(ylim=(1e-8, 2e-6)))

Ri_f histogram#

(
    np.log10(np.abs(baseline.Rif.cf.sel(Z=slice(-110, 0))))
    .reset_coords(drop=True)
    .hvplot.hist(bins=np.linspace(-2, 1, 100))
) * hv.VLine(np.log10(0.2)).opts(line_color="k")

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:771: RuntimeWarning: divide by zero encountered in log10
  result_data = func(*input_data)

ε-Jb histogram#

from xhistogram.xarray import histogram

jpdf = histogram(
    np.log10(baseline.eps),
    np.log10(0.2 * baseline.eps_Jb),
    density=True,
    bins=np.linspace(-10, -4, 100),
).compute()

jpdf.plot.contourf(levels=52, robust=True, cmap=mpl.cm.nipy_spectral)
dcpy.plots.line45()
_images/d71a0f78595ea225a4feff93c027ecb3b5375406401ada2cdfb1a610fadd0fb3.png 
from mom6_tools import wright_eos

mom6tao = baseline
mom6tao["α"] = wright_eos.alpha_wright_eos(mom6tao.thetao, mom6tao.so, p=0)
mom6tao["α"].attrs = {
    "standard_name": "sea_water_thermal_expansion_coefficient",
    "units": "kg m-3 C-1",
}
mom6tao["β"] = wright_eos.beta_wright_eos(mom6tao.thetao, mom6tao.so, p=0)
mom6tao["β"].attrs = {
    "standard_name": "sea_water_haline_contraction_coefficient",
    "units": "kg m-3",
}

(
    baseline.Kd_heat * baseline.α * baseline.Tz
    + baseline.Kd_heat * baseline.β * baseline.Sz
)

<xarray.DataArray (time: 507264, zi: 27, zl: 27)>
dask.array<chunksize=(8760, 26, 27), meta=np.ndarray>
Coordinates: (12/14)
  * time             (time) datetime64[ns] 1958-01-06T00:30:00 ... 2016-12-31...
    xh               float64 -140.0
    yh               float64 0.0625
    yq               float64 -0.0625
  * zi               (zi) float64 -230.8 -212.0 -194.4 -177.8 ... -5.0 -2.5 -0.0
    eucmax           (time) float64 dask.array<chunksize=(8760,), meta=np.ndarray>
    ...               ...
    en_mask          (time) bool False False False False ... False False False
    ln_mask          (time) bool False False False False ... False False False
    warm_mask        (time) bool True True True True ... True True True True
    cool_mask        (time) bool False False False False ... False False False
    enso_transition  (time) <U12 '____________' ... '____________'
  * zl               (zl) float64 -240.8 -221.4 -203.2 ... -6.25 -3.75 -1.25
Attributes:
    cell_measures:  area: areacello
    cell_methods:   area:mean zi:point yh:mean xh:mean time: mean
    long_name:      Total diapycnal diffusivity for heat at interfaces
    time_avg_info:  average_T1,average_T2,average_DT
    units:          m2 s-1

Mean SST, Jq#

def flatten(tree):
    varnames = next(iter(tree.children.values())).data_vars

    out = xr.Dataset()
    coord = xr.DataArray(list(tree.keys()), dims="node")
    for var in varnames:
        variables = (node[var] for name, node in tree.children.items())
        out[var] = xr.concat(variables, dim=coord)
    return out


sub = tree.map_over_subtree(
    lambda node: (
        node.cf[["Jq", "sea_surface_temperature", "ocean_vertical_heat_diffusivity"]]
        .pipe(np.abs)
        .cf.rename(
            {"sea_surface_temperature": "sst", "ocean_vertical_heat_diffusivity": "KT"}
        )
        .cf.sel(Z=slice(-60, -20))
        # .cf.mean("Z")
        .cf.groupby("time.month")
        .mean(["Z", "time"])
        .load()
    )
).pipe(flatten)

h = (
    sub.Jq.hvplot.line(by="node", legend=True)
    + sub.sst.hvplot.line(by="node", logy=True, legend=True)
).cols(1)
h

vertical = tree.map_over_subtree(
    lambda node: (
        node.cf[["Jq", "ocean_vertical_heat_diffusivity"]]
        .cf.rename({"ocean_vertical_heat_diffusivity": "KT"})
        .cf.sel(Z=slice(-120, -20))
        .groupby("time.season")
        .mean()
        .sel(season=["DJF", "MAM", "JJA", "SON"])
        .load()
    )
)

seasonal = {}
seasonal["Jq"] = mixpods.map_hvplot(
    lambda ds, name, muted: (
        ds["Jq"]
        .cf.rename({"Z": "depth"})
        .hvplot.line(y="depth", col="season", label=name, muted=muted)
    ),
    vertical.children,
)
seasonal["KT"] = mixpods.map_hvplot(
    lambda ds, name, muted: (
        ds["KT"]
        .cf.rename({"Z": "depth"})
        .hvplot.line(logx=True, y="depth", col="season", label=name, muted=muted)
    ),
    vertical.children,
)

def hvplot_facet(tree, varname, col, **kwargs):
    by = next(iter(tree.children.values()))[col].data
    handles = []
    for by_ in by:
        handles.append(
            hv.Overlay(
                [
                    node.ds.cf[varname]
                    .sel({col: by_})
                    .hvplot.line(label=name, **kwargs)
                    for name, node in tree.children.items()
                ]
            )
        )
    return hv.Layout(handles)


proc = vertical.map_over_subtree(lambda node: node.cf.rename({"Z": "Z"}))
layout = hv.Layout(
    [
        hvplot_facet(
            proc,
            varname="Jq",
            col="season",
        ),
        hvplot_facet(
            proc,
            varname="ocean_vertical_heat_diffusivity",
            col="season",
            logx=True,
            ylim=(1e-6, 1e-1),
        ),
    ]
)
layout.opts(
    hv.opts.Curve(invert_axes=True, xrotation=20, frame_width=160, frame_height=300),
    hv.opts.Overlay(legend_position="bottom_left", shared_axes=True, labelled=["x"]),
)

mixpods.map_hvplot(
    lambda ds, name, muted: (ds["Tz"] < 0)
    .sum("time")
    .hvplot.line(label=name, muted=muted, invert=True),
    tree,
)

np.log10(tree["kpp.lmd.004"]["Kd_heat"].resample(time="M").mean()).hvplot.quadmesh(
    y="zi"
)

mixpods.map_hvplot(
    lambda ds, name, muted: ds.ds.cf["ocean_vertical_heat_diffusivity"]
    .mean("time")
    .hvplot.line(label=name, muted=muted),
    tree.children,
).collate().opts(
    ylim=(1e-8, 1e1),
    legend_position="right",
    logx=True,
    invert_axes=True,
    frame_width=300,
    aspect=1 / 3,
)

WARNING:param.OverlayPlot02012: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.OverlayPlot02012: aspect value was ignored because absolute width and height values were provided. Either supply explicit frame_width and frame_height to achieve desired aspect OR supply a combination of width or height and an aspect value.

WARNING:param.OverlayPlot02072: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.OverlayPlot02072: aspect value was ignored because absolute width and height values were provided. Either supply explicit frame_width and frame_height to achieve desired aspect OR supply a combination of width or height and an aspect value.

mixpods.map_hvplot(
    lambda ds, name, muted: ds["Jq"].mean("time").hvplot.line(label=name, muted=muted),
    tree.children,
).collate().opts(
    legend_position="right", invert_axes=True, frame_width=300, aspect=1 / 3
)

Histograms: Ri, ε#

  • take out MLD?

handles = [
    mixpods.plot_distributions(tree, "chi", bins=np.linspace(-11, -4, 101), log=True),
    mixpods.plot_distributions(tree, "eps", bins=np.linspace(-11, -4, 101), log=True),
    mixpods.plot_distributions(
        tree, "ocean_vertical_heat_diffusivity", bins=np.linspace(-8, -1, 101), log=True
    ),
    # plot_distributions(tree, "Jq", bins=np.linspace(-1000, 0, 51), log=False),
    mixpods.plot_distributions(tree, "Rig_T", np.linspace(-0.5, 1.5, 61))
    * hv.VLine(0.25).opts(line_color="k"),
]

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:760: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:760: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)

hv.Layout(handles).opts("Overlay", frame_width=600).cols(2)

for node in ["baseline", "kpp.lmd.002", "kpp.lmd.003", "kpp.lmd.004"]:
    plt.figure()
    tree[node].sel(time=slice("2008-10-25", "2008-11-15")).ds.cf["ocean_vertical_heat_diffusivity"]
    .cf.sel(Z=slice(-120, None)).reset_coords(drop=True).cf.plot(
        cmap=mpl.cm.Spectral_r, size=3, aspect=5, norm=mpl.colors.LogNorm(1e-7, 1e-2)
    )

<Figure size 896x672 with 0 Axes>
_images/a53d29367e5dda09a0e1f6ac49d49c85e2f3fc6059252d26b38c3418d2c11c3e.png 
<Figure size 896x672 with 0 Axes>
_images/188bf451e3f81e0d5d66bc10b7cfbbcda8e63b0f7a13ae45b6ebd4a929c6e981.png 
<Figure size 896x672 with 0 Axes>
_images/b05865d1a8fbbb6c1b110ea191266f122a11440d655dcfe128666762ca0eb855.png 
<Figure size 896x672 with 0 Axes>
_images/dc2193a1caf3834b8809236cd1e5cdb814ef1548b3602f1a6feed6cb455942ee.png 
for node in ["baseline", "kpp.lmd.002", "kpp.lmd.003", "kpp.lmd.004"]:
    plt.figure()
    tree[node].sel(time=slice("2008-10-25", "2008-11-15"))["eps"].cf.sel(
        Z=slice(-120, None)
    ).reset_coords(drop=True).cf.plot(
        cmap=mpl.cm.Spectral_r, size=3, aspect=5, norm=mpl.colors.LogNorm(1e-9, 1e-5)
    )

<Figure size 896x672 with 0 Axes>
_images/afc6636caabd6f8f251314667803b93d38d79a6d71c5db2b1eae121ea814eb0e.png 
<Figure size 896x672 with 0 Axes>
_images/dab8621b127ca67dd6e59789fa32dc7a2503a41f8ec831816f553c78366f8144.png 
<Figure size 896x672 with 0 Axes>
_images/0924d5f87bf7108447438fec6abd2e7d33feb7cfd0c996b4af7b92af32a14f4a.png 
<Figure size 896x672 with 0 Axes>
_images/f42c5d1bd0bc2a1d39eafb97ca8a55a74a08c113156d791088e84e78effd7d70.png

Compare boundary layer depth#

hbl = (
    tree.dc.extract(["baseline", "kpp.lmd.004"])
    .dc.subset_nodes("KPP_OBLdepth")
    .dc.concatenate_nodes()
    .reset_coords(drop=True)
)

(
    hbl.groupby("time.hour").mean().hvplot.line(by="node", flip_yaxis=True)
    + hbl.groupby("time.hour").median().hvplot.line(by="node", flip_yaxis=True)
    + hbl.hvplot.hist(by="node", bins=np.arange(0, 50, 1), normed=1)
).cols(1)

Panel debug#

import panel as pn

print(layout)

:Layout
   .Overlay.I  :Overlay
      .Curve.TAO                             :Curve   [chi_bin]   (histogram_chi)
      .Curve.Baseline                        :Curve   [chi_bin]   (histogram_chi)
      .Curve.Epbl                            :Curve   [chi_bin]   (histogram_chi)
      .Curve.Kpp_full_stop_lmd_full_stop_002 :Curve   [chi_bin]   (histogram_chi)
      .Curve.Kpp_full_stop_lmd_full_stop_003 :Curve   [chi_bin]   (histogram_chi)
      .Curve.Kpp_full_stop_lmd_full_stop_004 :Curve   [chi_bin]   (histogram_chi)
   .Overlay.II :Overlay
      .Curve.TAO                             :Curve   [eps_bin]   (histogram_eps)
      .Curve.Baseline                        :Curve   [eps_bin]   (histogram_eps)
      .Curve.Epbl                            :Curve   [eps_bin]   (histogram_eps)
      .Curve.Kpp_full_stop_lmd_full_stop_002 :Curve   [eps_bin]   (histogram_eps)
      .Curve.Kpp_full_stop_lmd_full_stop_003 :Curve   [eps_bin]   (histogram_eps)
      .Curve.Kpp_full_stop_lmd_full_stop_004 :Curve   [eps_bin]   (histogram_eps)

def change_muted(value):
    present = layout.Overlay.I.Curve.children

    normalized = [n.replace(".", "_full_stop_") for n in value]
    [n for n in normalized if n in present]

    pass


def change_muted_test(value):
    present = layout.Overlay.I.Curve.children

    normalized = [n.replace(".", "_full_stop_") for n in value]
    [n for n in normalized if n in present]

    pass


names = list(tree.children.keys())
checkbox = pn.widgets.CheckBoxGroup(value=names, options=names, inline=True)
checkbox.jslink(layout, {"value": change_muted})

debug = pn.widgets.StaticText(name="debug", value="default")
checkbox.jslink(debug, dict(value=change_muted_test))

pn.Column(debug, checkbox, layout.cols(1))

bins = np.linspace(-0.5, 1.5, 61)
handles = [
    hvplot_step_hist(
        ds["Rig_T"].reset_coords(drop=True).cf.sel(Z=slice(-69, -29)),
        name=name,
        bins=bins,
    )
    for name, ds in datasets.items()
]
(hv.Overlay(handles) * hv.VLine(0.25).opts(line_color="k")).opts(ylabel="PDF")

ε vs Ri#

Possible reasons for difference

  1. My MLD seems deeper even though I use the same ΔT=0.15 threshold. It could be that they’ve used Akima splines to interpolate in the vertical

  2. ~I’ve filtered to DCL, so accounting for MLD and EUC movement. I’m not sure they did that.~

Only 𝜒pods between 29 and 69 m are used in this analysis as deeper 𝜒pods are more strongly influenced by the variability of zEUC than by surface forcing.

TODO

  • EUC strength is proportional to horizontal visc for 1° models

  • Add \(ε_χ\) for MOM6

  • Do for K

  • composite by TIW strength

  • start with 10min χpod data, then average to hourly.

  • make composites off the equator: look at strong off-equatorial du/dx; du/dz

f, ax = plt.subplots(
    2, len(datasets) // 2, sharex=True, sharey=True, constrained_layout=True
)
for axx, (name, ds) in zip(ax.flat, tree.children.items()):
    da = np.log10(ds["eps_ri"])
    # da["Rig_T_bins"] = pd.Index(
    #    [pd.Interval(np.log10(i.left), np.log10(i.right))
    #     for i in tao_gridded.Rig_T_bins.data]
    # )
    (
        da.sel(stat="mean")
        .sel(enso_transition_phase=["La-Nina cool", "El-Nino warm"])
        .plot(hue="enso_transition_phase", marker=".", ax=axx, add_legend=False)
    )
    axx.set_title(name)

    ticks = [0.04, 0.1, 0.25, 0.63, 1.6]
    # axx.set_yticks([-7.25, -7, -6.5, -6])
    axx.set_ylim([-7.5, -5.5])
    axx.set_xticks(np.log10(ticks))
    axx.set_xticklabels([f"{a:.2f}" for a in ticks])
    # axx.tick_params(axis='x', rotation=30)
ax[0, 0].set_ylabel("ε")
dcpy.plots.linex(np.log10(0.25), ax=ax.flat)
# dcpy.plots.clean_axes(ax)
f.set_size_inches((8, 4))
_images/c5c9326274e9719c4e1e803d16b84c2b6ea2a6fbdec8989d9a8caf9ed426ec83.png

Shear spectra#

Notes#

  • Now obvious difference between clockwise and counter-clockwise so I’ve chosen to just do the spectrum of the magnitude

  • think about tides

  • what is the “DATM coupling interval”

  • Kelvin’s KPP parameters

what are the momentum and heat fluxes over some resolvable timescale. And how do we compare to obs.

Frequency spectra @ Z= 50m#

DEPTHS = [-50, -60, -75]

f, ax = plt.subplots(
    len(DEPTHS), 1, sharex="col", sharey=True, squeeze=False, constrained_layout=True
)
for name in tree.children.keys():
    if name == "tropicheat":
        continue
    ds = tree[name].ds.reset_coords()
    u = ds.cf["sea_water_x_velocity"]
    v = ds.cf["sea_water_y_velocity"]
    S = (u + 1j * v).cf.sel(Z=DEPTHS, method="nearest").load()

    for axx, z in zip(ax, DEPTHS):
        dcpy.ts.PlotSpectrum(
            np.abs(S).cf.sel(Z=z, method="nearest"),
            multitaper=False,
            nsmooth=12,
            label=name,
            lw=0.75,
            ax=axx[0],
        )
        axx[-1].text(x=0.9, y=0.9, s=f"{z} m", transform=axx[-1].transAxes)
ax.ravel()[0].set_ylim([1e-9, None])
ax.ravel()[-1].legend()
dcpy.plots.clean_axes(ax)
f.set_size_inches((4, 6))
_images/2a9cc8145adfd44fee0e508cb0f5ea188843038db10e594c100c3a6b2d76b4a1.png 
f, ax = plt.subplots(
    len(DEPTHS), 1, sharex="col", sharey=True, squeeze=False, constrained_layout=True
)
for name in tree.children.keys():
    ds = tree[name].ds.reset_coords()
    u = ds.cf["sea_water_x_velocity"]
    v = ds.cf["sea_water_y_velocity"]
    S = (
        (u.cf.differentiate("Z") + 1j * v.cf.differentiate("Z"))
        .cf.sel(Z=DEPTHS, method="nearest")
        .load()
    )

    for axx, z in zip(ax, DEPTHS):
        dcpy.ts.PlotSpectrum(
            np.abs(S).cf.sel(Z=z, method="nearest"),
            multitaper=False,
            nsmooth=12,
            label=name,
            lw=0.75,
            ax=axx[0],
        )
        axx[-1].text(x=0.9, y=0.9, s=f"{z} m", transform=axx[-1].transAxes)
ax.ravel()[0].set_ylim([1e-9, None])
ax.ravel()[-1].legend()
dcpy.plots.clean_axes(ax)
f.set_size_inches((4, 6))
_images/a2bcd74b33b920df1c47a215f0e4fc97d6790adbfddaab01522a7c28e91d64aa.png

2D frequency-wavenumber spectra#

S = (u.cf.differentiate("Z") + 1j * v.cf.differentiate("Z")).sel(time="2010")
uniform = (
    S.cf.sel(Z=slice(-150, 0))
    .cf.interp(Z=np.arange(-150, -52, 5))
    .interpolate_na("time")
)
uniform

<xarray.DataArray (time: 8760, depth: 20)>
(-0.008227312937378883-0.014154188334941864j) ... (-0.030751610174775124+0.00...
Coordinates:
  * time     (time) datetime64[ns] 2010-01-01 ... 2010-12-31T23:00:00
  * depth    (depth) int64 -150 -145 -140 -135 -130 -125 ... -75 -70 -65 -60 -55
spec = xrft.power_spectrum(
    uniform,
    dim=(S.cf.axes["Z"][0], "time"),
    window="hann",
    window_correction=True,
)
sym = (
    2
    * 2
    * spec.isel(
        freq_depth=slice(uniform.cf.sizes["Z"] // 2 + 1, None),
        freq_time=slice(uniform.cf.sizes["time"] // 2 + 1, None),
    )
)

sym.mean("freq_depth").rolling(freq_time=5).mean(center=True, min_periods=5).plot(
    xscale="log", yscale="log"
)

[<matplotlib.lines.Line2D object at 0x2b8bababd5a0>]
_images/ecc7ae0a348800866276c7842528e70a44e323144d5be5b535dd3aaee85f8856.png 
import dcpy.ts

sym.plot(
    xscale="log",
    yscale="log",
    norm=mpl.colors.LogNorm(),
    robust=True,
    cmap=mpl.cm.turbo,
)
_images/b1fb54d7ba6fae2c436d1d4248449a4d9ad98d773a40a891796c349a396945f9.png

Vertical wavenumber spectra#

  • think about vertical migration of EUC

ds = tree["kpp.lmd.004"].ds

ds.S2.sel(time=ds.time.dt.month.isin([10, 11]))

<xarray.DataArray 'S2' (time: 36600, zi: 27)>
dask.array<chunksize=(1464, 26), meta=np.ndarray>
Coordinates: (12/14)
  * time             (time) datetime64[ns] 2003-10-01T00:30:00 ... 2027-11-30...
    xh               float64 -140.0
    yh               float64 0.0625
    yq               float64 -0.0625
  * zi               (zi) float64 -230.8 -212.0 -194.4 -177.8 ... -5.0 -2.5 -0.0
    eucmax           (time) float64 dask.array<chunksize=(1464,), meta=np.ndarray>
    ...               ...
    oni              (time) float32 0.2936 0.2936 0.2936 0.2936 ... nan nan nan
    en_mask          (time) bool False False False False ... False False False
    ln_mask          (time) bool False False False False ... False False False
    warm_mask        (time) bool True True True True ... True True True True
    cool_mask        (time) bool False False False False ... False False False
    enso_transition  (time) <U12 '____________' ... '____________'
Attributes:
    long_name:  $S^2$
    units:      s$^{-2}$
# f, ax = plt.subplots(1, 2, sharey=True, squeeze=False)
handles = []

for i, name in enumerate(tree.children.keys()):
    ds = tree[name].ds.reset_coords()
    u = ds.cf["sea_water_x_velocity"].load()
    v = ds.cf["sea_water_y_velocity"].load()
    S = (
        (u.cf.differentiate("Z") + 1j * v.cf.differentiate("Z"))
        .cf.sel(Z=slice(-200, 0))
        .cf.interp(Z=np.arange(-105, -52, 5))
    )
    Zname = S.cf.axes["Z"][0]
    spec = xrft.power_spectrum(
        S,
        dim=S.cf.axes["Z"],
        window="hann",
        window_correction=True,
    ).mean("time")
    sym = 2 * spec.isel({f"freq_{Zname}": slice(S.cf.sizes["Z"] // 2 + 1, None)})
    sym.name = "spectral density"
    handles.append(sym.hvplot.line(logx=True, logy=True, label=name))

hv.Overlay(handles).opts(legend_position="right")

# f, ax = plt.subplots(1, 2, sharey=True, squeeze=False)
handles = []

for i, name in enumerate(tree.children.keys()):
    ds = tree[name].ds.reset_coords()
    u = ds.cf["sea_water_x_velocity"].load()
    v = ds.cf["sea_water_y_velocity"].load()
    S = (
        (u.cf.differentiate("Z") + 1j * v.cf.differentiate("Z"))
        .cf.sel(Z=slice(-200, 0))
        .cf.interp(Z=np.arange(-105, -52, 5))
    )
    if "kpp" in name or "TAO" in name:
        S = S.sel(time=ds.time.dt.month.isin([10, 11]))
    Zname = S.cf.axes["Z"][0]
    spec = xrft.power_spectrum(
        S,
        dim=S.cf.axes["Z"],
        window="hann",
        window_correction=True,
    ).mean("time")
    sym = 2 * spec.isel({f"freq_{Zname}": slice(S.cf.sizes["Z"] // 2 + 1, None)})
    sym.name = "spectral density"
    handles.append(sym.hvplot.line(logx=True, logy=True, label=name))

hv.Overlay(handles).opts(legend_position="right")

f, ax = plt.subplots(1, 2, sharey=True, squeeze=False)
for name in tree.children.keys():
    ds = tree[name].ds.reset_coords()
    u = ds.cf["sea_water_x_velocity"].load()
    v = ds.cf["sea_water_y_velocity"].load()
    S = u.cf.differentiate("Z") + 1j * v.cf.differentiate("Z")
    dcpy.ts.PlotSpectrum(
        S.cf.sel(Z=-50, method="nearest").sel(time="2010"),
        multitaper=False,
        twoside=True,
        nsmooth=12,
        label=name,
        lw=0.75,
        ax=ax.ravel(),
    )
ax.ravel()[0].set_ylim([1e-10, None])
ax.ravel()[-1].legend()
f.set_size_inches((8, 4))
_images/fe8b29323164dac6185d3f4df5d5ab36ddabb6d0229fe307f54ccb71b72d9285.png

Seasonal mean heat flux#

remapped.Jq.load()

<xarray.DataArray 'Jq' (time: 174000, zeuc: 100)>
array([[        nan,         nan, -0.07841657, ...,         nan,
                nan,         nan],
       [        nan,         nan, -0.07973828, ...,         nan,
                nan,         nan],
       [        nan,         nan, -0.08094554, ...,         nan,
                nan,         nan],
       ...,
       [        nan, -0.07447129,         nan, ...,         nan,
                nan,         nan],
       [        nan, -0.07471326,         nan, ...,         nan,
                nan,         nan],
       [        nan, -0.07509062,         nan, ...,         nan,
                nan,         nan]])
Coordinates:
    latitude   float64 0.025
    longitude  float64 -140.0
  * time       (time) datetime64[ns] 1998-12-31T18:00:00 ... 2018-11-06T17:00:00
  * zeuc       (zeuc) float64 -247.5 -242.5 -237.5 -232.5 ... 237.5 242.5 247.5
(
    remapped.Jq.groupby("time.season")
    .mean()
    .reindex(season=["DJF", "MAM", "JJA", "SON"])
    .cf.plot.line(col="season")
)

<xarray.plot.facetgrid.FacetGrid>
_images/27e778991cb64aa8fa110aaf630c9fafb422676e51e9d35860ea1cc5ff1a5fe7.png

Below EUC mixing#

as a test for the background mixing params.

  • overflows

  • acc , gs

  • DWBC

  • ideal age metric?

  • vertically integrated heat budget term, global map

import os

import datatree

micro_zeuc = datatree.open_datatree(
    os.path.expanduser("~/datasets/microstructure/equix-tiwe-zeuc.average.nc")
).load()
for nodename, node in micro_zeuc.children.items():
    node["u"].attrs["standard_name"] = "sea_water_x_velocity"
    node["v"].attrs["standard_name"] = "sea_water_y_velocity"
    node["KT"].attrs["standard_name"] = "ocean_vertical_heat_diffusivity"
    node["ν"].attrs["standard_name"] = "ocean_vertical_momentum_diffusivity"

    del node["Rig_T"].attrs["standard_name"]
    del node["N2T"].attrs["standard_name"]
    # euc_mean[nodename] = node

%autoreload

newtree = mixpods.bin_to_euc_centered_coordinate(tree)

for nodename, _ in tree.children.items():
    tree[f"{nodename}/euc"] = newtree[f"{nodename}/euc"]

euc_mean = mixpods.average_euc(tree)
euc_mean.load()

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*
euc_mean.to_netcdf("mom6-euc-mean.nc")

for nodename, node in micro_zeuc.children.items():
    euc_mean[nodename] = node

h = {
    varname: mixpods.map_hvplot(
        lambda node, name, muted: (
            node.ds.cf[varname]
            .reset_coords(drop=True)
            .hvplot.line(
                ylabel=varname,
                label=name,
                logx=varname != "sea_water_x_velocity",
                invert=True,
            )
        ),
        euc_mean,
    )
    for varname in [
        "sea_water_x_velocity",
        "sea_water_potential_temperature",
        "chi",
        "eps",
        "ocean_vertical_heat_diffusivity",
        "ocean_vertical_momentum_diffusivity",
    ]
}
h["chi"].opts(ylim=(1e-9, 1e-4))
h["eps"].opts(ylim=(1e-9, 1e-4))
h["ocean_vertical_heat_diffusivity"].opts(ylim=(5e-7, 3))
h["ocean_vertical_momentum_diffusivity"].opts(ylim=(1e-6, 1e-1))

h2 = {
    varname: mixpods.map_hvplot(
        lambda node, name, muted: node.ds[varname]
        .reset_coords(drop=True)
        .hvplot.line(ylabel=varname, label=name, logx=False, invert=True),
        euc_mean,
    )
    for varname in ["Rig_T"]
}

hv.Layout(list(h.values()) + [h2["Rig_T"].opts(ylim=(0, 3))]).opts(
    hv.opts.Curve(frame_width=150, frame_height=300, xlim=(-150, 150)),
    hv.opts.Layout(shared_axes=True),
    hv.opts.Overlay(show_legend=True, show_grid=True, legend_position="right"),
    *mixpods.HV_TOOLS_OPTIONS,
    *mixpods.PRESENTATION_OPTS,
).cols(4)

Ri = (
    tree.dc.extract_leaf("euc")
    .dc.subset_nodes(["Rig_T"])
    .sel(time=slice("1992", "2017"))
    .load()
)

mixpods.plot_distributions(
    Ri.dc.reorder_nodes(["TAO", ...]),
    "Rig_T",
    bins=np.arange(0, 5, 0.25),
    log=False,
    subset=False,
).opts(hv.opts.Overlay(legend_position="right"))

# tree["new_baseline.hb/euc"]["Rig_T"].sel(zeuc=slice(-120, -30)).hvplot.hist(bins=[0, 0.2, 0.5, 0.8, 1, 2, 3, 4, 5])
# tree["baseline/euc"]["Rig_T"].sel(zeuc=slice(-120, -30)).hvplot.hist(bins=[0, 0.2, 0.5, 0.8, 1, 2, 3, 4, 5])

catalog.search(frequency="monthly").df
casename stream path baseline levels frequency variables shortname description
0 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [N2_int, agessc, h, rhopot0, so, thetao, uhGM,... baseline.001 baseline
1 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... hm /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [Heat_PmE, KE, Rd_dx, SSH, SSU, SSV, T_adx_2d,... baseline.001 baseline
2 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [Rd_dx, SSH, SSU, SSV, mass_wt, mlotst, oml, o... baseline.001 baseline
3 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [h, so, thetao, uhGM, uhml, umo, uo, vhGM, vhm... baseline.epbl.001 ePBL
4 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [SSH, SSU, SSV, mlotst, sos, speed, tos] baseline.epbl.001 ePBL
5 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... wci /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [T_advection_xy, T_lbdxy_cont_tendency, Th_ten... baseline.epbl.001 ePBL
6 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [h, so, thetao, uhGM, uhml, umo, uo, vhGM, vhm... baseline.kpp.lmd.002 KPP Ri0=0.5
7 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [SSH, SSU, SSV, mlotst, oml, sos, speed, tos] baseline.kpp.lmd.002 KPP Ri0=0.5
8 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... wci /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [KPP_NLT_temp_budget, T_advection_xy, T_lbdxy_... baseline.kpp.lmd.002 KPP Ri0=0.5
9 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [h, so, thetao, uhGM, uhml, umo, uo, vhGM, vhm... baseline.kpp.lmd.003 KPP Ri0=0.5, Ric=0.2,
10 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [SSH, SSU, SSV, mlotst, oml, sos, speed, tos] baseline.kpp.lmd.003 KPP Ri0=0.5, Ric=0.2,
11 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... wci /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [KPP_NLT_temp_budget, T_advection_xy, T_lbdxy_... baseline.kpp.lmd.003 KPP Ri0=0.5, Ric=0.2,
12 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [h, so, thetao, uhGM, uhml, umo, uo, vhGM, vhm... baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
13 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [SSH, SSU, SSV, mlotst, oml, sos, speed, tos] baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
14 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... wci /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [KPP_NLT_temp_budget, T_advection_xy, T_lbdxy_... baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
15 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [N2_int, agessc, h, rhopot0, so, thetao, uhGM,... new_baseline.hb KD=0, KV=0
16 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... hm /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [Heat_PmE, KE, KPP_NLT_temp_budget, Rd_dx, SSH... new_baseline.hb KD=0, KV=0
17 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [N2_int, agessc, h, rhopot0, so, thetao, uhGM,... new_baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
18 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... hm /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [Heat_PmE, KE, KPP_NLT_temp_budget, Rd_dx, SSH... new_baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
19 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [N2_int, agessc, h, rhopot0, so, thetao, uhGM,... new_baseline.kpp.lmd.005 KPP ν0=2.5, Ri0=0.5
20 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... hm /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [Heat_PmE, KE, KPP_NLT_temp_budget, Rd_dx, SSH... new_baseline.kpp.lmd.005 KPP ν0=2.5, Ri0=0.5 
monthly = catalog.search(frequency="monthly", variables="uo").to_datatree()

--> The keys in the returned dictionary of datasets are constructed as follows:
	'shortname/stream'
100.00% [8/8 00:00<00:00] 
uo = (
    monthly.drop_nodes(
        set(monthly)
        - set(
            tree.dc.rename_nodes(
                {"baseline": "baseline.001", "kpp.lmd.004": "baseline.kpp.lmd.004"}
            )
        )
    )
    .dc.extract_leaf("h")
    .sel(xq=-140, method="nearest")
    .sel(yh=slice(-10, 10))
    .sel(time=slice("0046-01-01", "0058-01-01"))
    # .mean("time")
).load()

(
    hv.Layout(
        [
            node["uo"]
            .where(node["uo"] < 10)
            .reset_coords(drop=True)
            .hvplot.quadmesh(x="yh", clim=(-0.01, 0.01), cmap="coolwarm", title=name)
            for name, node in uo.mean("time").children.items()
        ]
    )
    .cols(2)
    .opts(
        hv.opts.QuadMesh(invert_yaxis=True, ylim=(2000, 0)),
        *mixpods.PRESENTATION_OPTS,
        *mixpods.HV_TOOLS_OPTIONS,
    )
)

Test out available variables#

  1. KPP_Kheat is non-zero only in KPP_OBL

(mom6140.Tflx_dia_diff * 1025 * 4000).cf.plot(robust=True)

<matplotlib.collections.QuadMesh>
_images/380ae0927607e64c1f1e455f3c2bcd9a0a70306e4e3510a888f8b10c5fd4535d.png 
mom6140.Kv_u.cf.plot(norm=mpl.colors.LogNorm(1e-6, 1e-2))

<matplotlib.collections.QuadMesh>
_images/3b26c61324b5edb9b740b26809160c8c09dff78025129659a54429702b7bc67b.png 
mom6140.Kd_heat.cf.plot(norm=mpl.colors.LogNorm(1e-6, 1e-2))
mom6140.KPP_OBLdepth.plot(color="r")

plt.figure()
mom6140.KPP_kheat.cf.plot(norm=mpl.colors.LogNorm(1e-6, 1e-2))
mom6140.KPP_OBLdepth.plot(color="r")

plt.figure()
(mom6140.KPP_kheat - mom6140.Kd_heat).cf.plot(
    robust=True, cmap=mpl.cm.mpl.cm.Reds_r, vmax=0
)

<matplotlib.collections.QuadMesh>
_images/3e83df3f9fea04390152de895b1d57178bcf8b50cf8f974de09d4ccaf8badbd6.png _images/f0f31fc31fc6c5f88350878aaf6306750183ad59680e27b9c7a37e7efb9042ce.png _images/2daecb46b41b3924e5589d365513732343962be989936b3c21f36dddb2677df9.png

Experiment with manual combining#

from intake.source.utils import reverse_format

years = []
tiles = []
for file in files[:10]:
    p = pathlib.Path(file)
    params = reverse_format("__{year}.nc.{tile}", p.stem + p.suffix)
    years.append(params["year"])
    tiles.append(params["tile"])
years, tiles

import toolz as tlz

bytile = {}
for tile, v in tlz.groupby(tileid, files).items():
    bytile[tile] = xr.concat(read_raw_files(v, parallel=True), dim="time")


print("\n".join([hash(ds.yh.data.tolist()) for ds in list(bytile.values())]))

from functools import partial


def hash_coords(ds, axis):
    dims = ds.cf.axes[axis]
    data = np.concatenate([ds[dim].data for dim in dims])
    return hash(tuple(data))


grouped = bytile

for axis, concat_axis in [("X", "Y"), ("Y", "X")]:
    grouped = tlz.groupby(partial(hash_coords, axis=axis), grouped.values())
    grouped = {
        k: cfconcat(v, axis=concat_axis, join="exact") for k, v in grouped.items()
    }
    break
grouped

cfconcat(list(grouped.values()), "X")

combined = xr.combine_by_coords(list(bytile.values()), combine_attrs="override")


def raise_if_bad_index(combined):
    bad = []
    for idx in combined.indexes:
        index = combined.indexes[idx]
        if not index.is_monotonic or index.has_duplicates:
            bad.append(idx)
    if bad:
        raise ValueError(f"Indexes {idx} are either not monotonic or have duplicates")


def tileid(path):
    p = pathlib.Path(path)
    # print(p.suffix)
    # params = reverse_format("__{year}.nc.{tile}", p.stem + p.suffix)
    return int(p.suffix[1:])  # params["tile"]
    # years.append(params["year"])
    # tiles.append(params["tile"])


sorted_files = sorted(files, key=tileid)
for tile, files in groupby(sorted_files, tileid):
    print(tile, len(list(files)))

Test out ONI#

Phase labeling#

oni = pump.obs.process_oni().sel(time=slice("2005-Sep", None))
enso_transition = mixpods.make_enso_transition_mask(oni)
mixpods.plot_enso_transition(oni, enso_transition)

Replicate ONI calculation#

close enough!

ersst = xr.tutorial.open_dataset("ersstv5")

monthlyersst = (
    ersst.sst.cf.sortby("Y")
    .cf.sel(latitude=slice(-5, 5), longitude=slice(360 - 170, 360 - 120))
    .cf.mean(["X", "Y"])
    .resample(time="M")
    .mean()
    .load()
)

expected = mixpods.calc_oni(monthlyersst)

actual = pump.obs.process_oni()
actual.plot()
expected.plot()

(actual - expected).plot()

[<matplotlib.lines.Line2D>]
_images/77d06dddd6387c72a89ff81237d73e61d968f573f74e1176bd004cbc9a48d694.png

Comparing to Warner and Moum code#

chipod.eps.sel(time="2007-01-14 10:35:00", method="nearest").load()

<xarray.DataArray 'eps' (depth: 5)>
nan 3.929e-07 2.846e-07 2.234e-07 4.491e-07
Coordinates:
    time        datetime64[ns] 2007-01-14T11:00:00
  * depth       (depth) float64 -69.0 -59.0 -49.0 -39.0 -29.0
    timeSeries  (depth) float64 69.0 59.0 49.0 39.0 29.0
    lat         (depth) float64 0.0 0.0 0.0 0.0 0.0
    lon         (depth) float64 -140.0 -140.0 -140.0 -140.0 -140.0
Attributes:
    long_name:              turbulence dissipation rate
    standard_name:          specific_turbulent_kinetic_energy_dissipation_in_...
    ncei_name:              turbulent kinetic energy (TKE) dissipation rate
    units:                  W kg-1
    FillValue:              -9999
    valid_min:              1e-12
    valid_max:              9.999999999999999e-06
    coverage_content_type:  physicalMeasurement
    grid_mapping:           crs
    source:                 inferred from fast thermistor spectral scaling
    references:             Moum J.N. and J.D. Nash, Mixing measurements on a...
    cell_methods:           depth: point, time: mean
    platform:               mooring
    instrument:             chipod
tao_gridded.enso_transition.loc[{"time": slice("2015-12", "2016-01")}]

<xarray.DataArray 'enso_transition' (time: 1488)>
'El-Nino warm' 'El-Nino warm' 'El-Nino warm' ... 'El-Nino warm' 'El-Nino warm'
Coordinates: (12/13)
    deepest             (time) float64 -300.0 -300.0 -300.0 ... -300.0 -300.0
    eucmax              (time) float64 -135.0 -140.0 -135.0 ... -140.0 -135.0
    latitude            float32 0.0
    longitude           float32 -140.0
    mld                 (time) float64 -9.0 -8.0 -8.0 -8.0 ... nan nan nan nan
    mldT                (time) float64 -10.0 -8.0 -8.0 -8.0 ... nan nan nan nan
    ...                  ...
    shallowest          (time) float64 -5.0 -5.0 -5.0 -5.0 ... -10.0 -10.0 -10.0
  * time                (time) datetime64[ns] 2015-12-01 ... 2016-01-31T23:00:00
    oni                 (time) float32 2.53 2.53 2.53 2.53 ... 2.53 2.53 2.53
    warm_mask           (time) bool True True True True ... False False False
    cool_mask           (time) bool False False False False ... True True True
    enso_transition     (time) <U12 'El-Nino warm' ... 'El-Nino warm'
Attributes:
    description:  Warner & Moum (2019) ENSO transition phase; El-Nino = ONI >...

Checking background visc, diff#

staticfile = "/glade/campaign/cgd/oce/projects/pump/cesm//gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb/run/*static*.nc"
static = xr.open_mfdataset(staticfile).isel(time=1).squeeze()

static.Kd_bkgnd.cf.sel(Z=200, method="nearest").hvplot.quadmesh(
    cmap="spectral_r", cnorm="log", clim=(5e-7, 2e-4)
)

ds = static.cf.sel(xh=5, yh=-25, method="nearest")

(
    ds.Kd_bkgnd.load().hvplot.line(label="Kd_bkgnd")
    * ds.Kv_bkgnd.load().hvplot.line(label="Kv_bkgnd")
).opts(hv.opts.Curve(logx=True, ylim=(1e-7, 1e-4), invert_axes=True)).opts(
    legend_position="right", show_legend=True, width=400, height=400
)

ds = static.cf.sel(xh=-140, yh=0, method="nearest")

(
    ds.Kd_bkgnd.load().hvplot.line(label="Kd_bkgnd")
    * ds.Kv_bkgnd.load().hvplot.line(label="Kv_bkgnd")
).opts(hv.opts.Curve(logx=True, ylim=(1e-7, 1e-4), invert_axes=True)).opts(
    legend_position="right", show_legend=True, width=400, height=400
)

opottempdiff vs Tflx_dia_diff#

ds = mixpods.load_mom6_sections(
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods"
)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yq', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yq', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(

grid = xgcm.Grid(ds, coords={("Z",): {"left": "zi", "center": "zl"}})
grid

<xgcm.Grid>
('Z',) Axis (periodic, boundary=None):
  * left     zi --> center
  * center   zl --> left

dz = ds.zi.diff("zi").rename({"zi": "zl"}).variable
dz

<xarray.Variable (zl: 36)>
array([42.79, 38.5 , 34.87, 31.76, 29.1 , 26.79, 24.77, 23.  , 21.42,
       20.02, 18.76, 17.61, 16.56, 15.59, 14.69, 13.85, 13.06, 12.29,
       11.54, 10.81, 10.08,  9.37,  8.66,  7.97,  7.28,  6.61,  5.95,
        5.29,  4.65,  4.01,  3.38,  2.77,  2.5 ,  2.5 ,  2.5 ,  2.5 ])
Attributes:
    cartesian_axis:  Z
    long_name:       Interface pseudo-depth, -z*
    positive:        up
    units:           meter
ds.Tflx_dia_diff.load()
ds.opottempdiff.load()

<xarray.DataArray 'opottempdiff' (time: 218424, zl: 37)>
0.004551 0.006139 0.01125 0.01764 0.01268 ... 39.39 62.65 12.85 10.76 -125.7
Coordinates: (12/13)
  * time             (time) datetime64[ns] 2003-01-07T00:30:00 ... 2028-01-06...
    xh               float64 -140.0
    yh               float64 0.0625
    yq               float64 -0.0625
  * zl               (zl) float64 -547.8 -502.4 -461.8 ... -6.25 -3.75 -1.25
    eucmax           (time) float64 -114.5 -114.5 -114.5 ... -114.5 -114.5
    ...               ...
    oni              (time) float32 nan nan nan nan nan ... nan nan nan nan nan
    en_mask          (time) bool False False False False ... False False False
    ln_mask          (time) bool False False False False ... False False False
    warm_mask        (time) bool True True True True ... True True True True
    cool_mask        (time) bool False False False False ... False False False
    enso_transition  (time) <U12 '____________' ... '____________'
Attributes:
    cell_measures:  volume: volcello area: areacello
    cell_methods:   area:mean zl:sum yh:mean xh:mean time: mean
    long_name:      Tendency of sea water potential temperature expressed as ...
    standard_name:  tendency_of_sea_water_potential_temperature_expressed_as_...
    time_avg_info:  average_T1,average_T2,average_DT
    units:          W m-2
integ = (ds.opottempdiff.isel(zl=slice(1, None)) * dz).cumsum("zl")

(1028 * 4200 * ds["Tflx_dia_diff"]).isel(time=0).differentiate("zi").plot()
ds["opottempdiff"].isel(time=0).plot()
# integ.isel(time=0).plot()

[<matplotlib.lines.Line2D object at 0x2ac9c184fbb0>]
_images/73b000a016a621f96408dc93538661e016d07bcd6a96cc2fa316694bbdcfdc93.png

Debugging total visc, diff:#

Bug report: NCAR/MOM6#238

Kv_u should it match the sum of in diffusivities

(
    ds.Kd_heat.load().mean("time").hvplot.line(label="Kd_heat")
    * ds.Kv_u.load().mean("time").hvplot.line(label="Kv_u")
    * ds.Kv_v.load().mean("time").hvplot.line(label="Kv_v")
).opts(hv.opts.Curve(logx=True, ylim=(1e-7, None), invert_axes=True)).opts(
    legend_position="right", show_legend=True, width=400, height=400
)

Check Kd, Kv contributions#

What is happening with viscosity at 80-m? Lateral averaging? vertical filter?

  • somehow laundered through the tri-diagonal solver

need internal tidal dissipation; CVMix_Tidal; This is near-bottom

Need Kv_itides

moor = xr.open_dataset(
    "/glade/u/home/dcherian/scratch/"
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods/run/"
    # "kd0kv0/"
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods.mom6.TAO_140W__0046.nc.0460",
    use_cftime=True,
).pipe(preprocess_mom6_sections)

mixpods.decompose_diff_visc(moor)

returning average_T1
returning average_T2
returning average_DT
returning time_bnds
returning average_T1
returning average_T2
returning average_DT
returning time_bnds

moor = xr.open_dataset(
    "/glade/u/home/dcherian/scratch/"
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods/run/"
    "kd0kv0/"
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods.mom6.TAO_140W__0046.nc.0460",
    use_cftime=True,
)
mixpods.decompose_diff_visc(moor)

smoothing of shear only?

moor = xr.open_dataset(
    "/glade/u/home/dcherian/scratch/"
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods/run/"
    "kd0kv0/"
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods.mom6.TAO_140W__0046.nc.0493",
    use_cftime=True,
)
mixpods.decompose_diff_visc(moor)

KD=0, KV=0 (kd0kv0/)#

moor = xr.open_dataset(
    "/glade/u/home/dcherian/scratch/"
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods/run/"
    "kd0kv0/"
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods.mom6.TAO_140W__0046.nc.0428",
    use_cftime=True,
)
mixpods.decompose_diff_visc(moor)

KD=2e-5, KV=1e-4 (kddefkvdef)#

moor = xr.open_dataset(
    "/glade/u/home/dcherian/scratch/"
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods/run/"
    "kddefkvdef/"
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods.mom6.TAO_140W__0046.nc.0428",
    use_cftime=True,
)
mixpods.decompose_diff_visc(moor)

Extra code#